function dy = fieldeqns(t,y,epsilont, epsilon, chi, kappa1, kappa2, Delta_a, Delta_b, a, b, c)


dy = zeros(2,1);

chi = -chi*((-1)^a + (-1)^b + (-1)^c);
eps = interp1(epsilont,epsilon,t); % Interpolate the data set (gt,g) at time t

dy(1) = -sqrt(kappa1)*eps - 1j*(Delta_a + chi)*y(1) - 0.5*kappa1*y(1) - 0.5*sqrt(kappa1*kappa2)*y(2);
dy(2) = -sqrt(kappa2)*eps - 1j*(Delta_b + chi)*y(2) - 0.5*kappa2*y(2) - 0.5*sqrt(kappa1*kappa2)*y(1);

%% solution corresponding to perfect parity measurement

% chi = -chi*((-1)^a + (-1)^b + (-1)^c);
% 
% if chi == 1
%     chi = -3;
% end
% 
% if chi == -1
%     chi = 3;
% end
% 
% eps = interp1(epsilont,epsilon,t); % Interpolate the data set (gt,g) at time t
% 
% dy(1) = -sqrt(kappa1)*eps - 1j*(Delta_a + chi)*y(1) - 0.5*kappa1*y(1) - 0.5*sqrt(kappa1*kappa2)*y(2);
% dy(2) = -sqrt(kappa2)*eps - 1j*(Delta_b + chi)*y(2) - 0.5*kappa2*y(2) - 0.5*sqrt(kappa1*kappa2)*y(1);

end